*	See Gibbs estimation of microstructure models: Teaching notes
	November, 2006
	Joel Hasbrouck
____________________________________________________________________________________________________

	RollGibbs: macro to estimate Roll model with Gibbs sampler

	RollGibbs(
		dsIn, 				input dataset with price variable 'p'
		qOut=qOut, 			output dataset containing simulated q's
		parmOut=parmOut, 	output dataset containing simulated parameters
		nSweeps=100,		number of sweeps
		qDraw=1, 			Draw (simulate) the q's (set qDraw=0 otherwise)
		varuDraw=1, 		Draw varu
		regDraw=1, 			Draw the regression coefficient (c)
		varuStart=0.001, 	Starting value for varu
		cStart=.01, 		Starting value for c
		printLevel=0)

____________________________________________________________________________________________________;


%let Infinity=1e30;
%let eps=1e-30;

%macro RollGibbs(dsIn, qOut=qOut, parmOut=parmOut, nSweeps=100,
	qDraw=1, varuDraw=1, regDraw=1, varuStart=0.001, cStart=.01, printLevel=0,
	cLower=0, cUpper=&Infinity);
proc iml;
	start main;
	reset printadv=1;
	call streaminit(1234);	*	Initialize the random number generators;
	reset storage=this.imlstor; load;	*	Reload necessary subroutines;
	show modules;

	*	Read in data;
	use &dsIn;
	if &qDraw=1 then do;
		read all var {p};
		q = j(nrow(p),1,.);
	end;
	else read all var {p q};
	nObs = nrow(p);

	*	Create output datasets for simulated qs and parameters;
	create &qOut var {sweep t q};
	create &parmOut var {sweep sdu c};

	dp = p[2:nObs] - p[1:(nObs-1)];

	if &qDraw then do; 	*	Initialize qs to price sign changes;
		qInitial = {1} // sign(dp);
		qInitial = qInitial # (q^=0);	*	Only initialize nonzero elements of q;
		q = qInitial;
	end;

	varu = &varuStart;
	c = &cStart;

	do sweep=1 to &nSweeps;

		if mod(sweep,1000)=0 then print sweep [label='' rowname='Sweep:'];

		dq =  q[2:nObs] - q[1:(nObs-1)];

		if &regDraw then do;
			priorMu = j(1,1,0);
			postMu = priorMu;
			priorCov = j(1,1,1);
			postCov = priorCov;
			rc = BayesRegressionUpdate(priorMu, priorCov, dp, dq, varu, postMu, postCov);
			if &printLevel>=2 then print postMu postCov;
			c = mvnrndT(postMu, postCov, &cLower, &cUpper);	*	Draw c (nonnegative);
			if &printLevel>=2 then print c;
		end;

		if &varuDraw then do;
			u = dp - c*dq;
			priorAlpha = 1.e-12;
			priorBeta = 1.e-12;
			postAlpha = .;
			postBeta = .;
			rc = BayesVarianceUpdate(priorAlpha, priorBeta, u, postAlpha, postBeta);
			x = (1/postBeta) * rand('gamma',postAlpha);
			varu = 1/x;
			sdu = sqrt(varu);
			if &printLevel>=2 then print varu;
		end;

		if &qDraw then do;
			qDrawPrintLevel = 0;
			call qDraw(p, q, c, varu, qDrawPrintLevel);
/*			call qDrawSlow(p, q, c, varu, qDrawPrintLevel);*/
			setout &qOut;
			mOut = j(nObs,1,sweep) || colvec(1:nObs) || colvec(q);
			append from mOut;
		end;

		if &regDraw | &varuDraw then do;
			setout &parmOut;
			append var _all_;
		end;
	end;
	finish main;

run;
quit;
%mend RollGibbs;
